{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Open Source GIS Pre-Processing Tutorial\n",
    "\n",
    "This notebook will rely on the WRF-Hydro GIS Pre-processing tools, found here:  \n",
    "* https://github.com/NCAR/wrf_hydro_gis_preprocessor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "1. [Create domain boundary shapefile](#1.-Create-domain-boundary-shapefile)<br>\n",
    "2. [Build GeoTiff raster from a WPS Geogrid file](#2.-Build-GeoTiff-raster-from-a-WPS-Geogrid-file)<br>\n",
    "3. [Building the hydrologic routing grids, aka the \"routing stack\"](#3.-Building-the-hydrologic-routing-grids,-aka-the-\"routing-stack\")<br>\n",
    "4. [Understanding the outputs](#4.-Understanding-the-outputs)\n",
    "5. [Examine outputs of GIS pre-processor](#5.-Examine-outputs-of-GIS-pre-processor)\n",
    "6. [Optional - Build Non-NWM WRF-Hydro Configurations of the Croton, NY Test Case](#6.-[Optional]-Build-Non-NWM-WRF-Hydro-Configurations-of-the-Croton,-NY-Test-Case)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### First, set the file paths for inputs and outputs\n",
    "\n",
    "Throughout this exercise, we will use python variables to store directory paths and other variables. However, we will call all GIS Pre-processing functionality as though it were on the command line. This is done by adding `!` syntax before the command-line syntax, to execute the line using bash.\n",
    "\n",
    "In this cell, Python variables are created that point to the file paths of the test-case data and an output directory is defined to store the data created by these tools."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Import python core modules\n",
    "import os\n",
    "import shutil\n",
    "\n",
    "# Set root directory for GIS lesson\n",
    "gis_data_folder = \"/home/docker/GIS_Training\"\n",
    "\n",
    "# Change the directory to the GIS_Training directory and get current working directory\n",
    "os.chdir(gis_data_folder)\n",
    "cwd = os.getcwd()\n",
    "\n",
    "# Set paths to known input and output directories and files\n",
    "data_folder = os.path.join(cwd, 'Croton_Lambert')\n",
    "in_geogrid = os.path.join(data_folder, 'geo_em.d01.nc')\n",
    "output_folder = os.path.join(cwd, 'Outputs')\n",
    "\n",
    "# Clear any outputs from previous runs by deleting (if necessary) and re-creating the output directory\n",
    "if os.path.exists(output_folder):\n",
    "    shutil.rmtree(output_folder)\n",
    "os.mkdir(output_folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Create domain boundary shapefile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The tool `Create_Domain_Boundary_Shapefile.py` takes a WRF (WPS) output file, aka \"Geogrid file\", and creates a polygon shapefile that defines the boundary of the domain as a single rectangular polygon in projected coordinates. The script will read metadata in the geogrid file and the output shapefile will be in the projection of the WRF domain. The unstaggered grid, or \"Mass\" grid (e.g. \"HGT_M\" variable), is used as the routing grid domain by WRF-Hydro."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Request help message from the script\n",
    "This is an example of the syntax for calling the `Create_Domain_Boundary_Shapefile.py` tool on the command line. By following the tool with `-h` or `--help`, we are able to call the help arguement which explains the purpose of the tool, shows the different arguements we can use for this tool, as well as the descriptions for each arguement. \n",
    "\n",
    "When the tool is run in the terminal, it is not necessary to use the exclamation point. In Jupyter Notebook, we can execute command-line syntax using \"!\". "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Script initiated at Tue Nov  3 00:15:03 2020\n",
      "usage: Create_Domain_Boundary_Shapefile.py [-h] -i IN_NC -o OUT_DIR\n",
      "\n",
      "This tool takes an WRF Geogrid file and creates a single polygon shapefile\n",
      "that makes up the boundary of the domain of the M-grid (HGT_M, for example).\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help  show this help message and exit\n",
      "  -i IN_NC    Path to WPS geogrid (geo_em.d0*.nc) file or WRF-Hydro\n",
      "              Fulldom_hires.nc file.\n",
      "  -o OUT_DIR  Output directory.\n"
     ]
    }
   ],
   "source": [
    "# Execute script on the command-line, requesting tool help (parameter -h)\n",
    "! python Create_Domain_Boundary_Shapefile.py -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "You will see in the messages above that the tool provides a brief explanation of the expected input and output parameters. This tool requires a geogrid file as input (`-i`) and a directory to write the outputs into (`-o`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "#### Execute the script using command-line syntax\n",
    "Now that we know what arguements are needed for this tool, we can enter those arguements and run the tool. For this script, we only need to specify the file path to the WPS geogrid file and an output folder to save the result. The result of this tool is a shapefile that shows the geographic boundary of the domain, as defined in the geogrid file. \n",
    "\n",
    "When running this tool in Jupyter, we can use brackets around our python variable names, and Jupyter will substitute the variable values when executing the syntax. This is akin to using an environment variable on the command-line. For the sake of repeatability, we also print the full syntax for reference. This can be copied into the terminal if desired."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Command to run:\n",
      "\n",
      "python Create_Domain_Boundary_Shapefile.py \\\n",
      "\t -i /home/docker/GIS_Training/Croton_Lambert/geo_em.d01.nc \\\n",
      "\t -o /home/docker/GIS_Training/Outputs\n",
      "\n",
      "Script initiated at Tue Nov  3 00:15:04 2020\n",
      "WPS netCDF projection identification initiated...\n",
      "    Map Projection: Lambert Conformal Conic\n",
      "    Using MOAD_CEN_LAT for latitude of origin.\n",
      "    Using Standard Parallel 2 in Lambert Conformal Conic map projection.\n",
      "    Geo-referencing step completed without error in  0.05 seconds.\n",
      "    Created projection definition from input NetCDF GEOGRID file.\n",
      "      ESRI Shapefile driver is available.\n",
      "  Done producing output vector polygon shapefile in  0.00 seconds\n",
      "  Output shapefile: /home/docker/GIS_Training/Outputs/geo_em.d01_boundary.shp\n",
      "Process completed in 0.06 seconds.\n"
     ]
    }
   ],
   "source": [
    "# Print information to screen for reference\n",
    "print('Command to run:\\n')\n",
    "print('python Create_Domain_Boundary_Shapefile.py \\\\\\n\\t -i {0} \\\\\\n\\t -o {1}\\n'.format(in_geogrid, output_folder))\n",
    "\n",
    "# Run the script with required parameters\n",
    "! python Create_Domain_Boundary_Shapefile.py -i {in_geogrid} -o {output_folder}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "----\n",
    "The messages returned by the tool can be quite useful. You will see the coordinate system information printed to the screen and any other progress messages. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Build GeoTiff raster from a WPS Geogrid file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The tool `Build_GeoTiff_From_Geogrid_File.py` is a program to export variables from a WRF-Hydro input file (geogrid or Fulldom_hires) file to an output raster format, with all spatial and coordinate system metadata. If a 3-dimensional variable is selected, individual raster bands will be created in the output raster for each index in the 3rd dimension. If a 4-dimensional variable is selected, the first index in the 4th dimension will be selected and the variable will be treated as a 3-dimensional variable described above.\n",
    "\n",
    "This tool is handy for performing a quick vizualization using GIS or othe software to examine the contents of the WRF-Hydro input file and overlay these grids with other goespatial data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tool takes three input parameters: an input Geogrid or Fulldom_hires netCDF file (`-i`), a variable name (`-v`), and an output GeoTiff raster file (`-o`) that the tool will create. For this example, we will export the variable \"HGT_M\", or surface elevation in meters above sea level.\n",
    "\n",
    "#### Request help message from the script"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Script initiated at Tue Nov  3 00:15:06 2020\n",
      "usage: Build_GeoTiff_From_Geogrid_File.py [-h] -i IN_NC [-v VARIABLE]\n",
      "                                          [-o OUT_FILE]\n",
      "\n",
      "This is a program to export >=2D variables from a WRF-Hydro input file\n",
      "(geogrid or Fulldom_hires) file to an output raster format, with all spatial\n",
      "and coordinate system metadata. If a 3-dimensional variable is selected,\n",
      "individual raster bands will be created in the output raster for each index in\n",
      "the 3rd dimension. If a 4-dimensional variable is selected, the first index in\n",
      "the 4th dimension will be selected and the variable will be treated as a\n",
      "3-dimensional variable described above.\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help   show this help message and exit\n",
      "  -i IN_NC     Path to WPS geogrid (geo_em.d0*.nc) file or WRF-Hydro\n",
      "               Fulldom_hires.nc file.\n",
      "  -v VARIABLE  Name of the variable in the input netCDF file. default=HGT_M\n",
      "  -o OUT_FILE  Output GeoTiff raster file.\n"
     ]
    }
   ],
   "source": [
    "# Get script help information\n",
    "! python Build_GeoTiff_From_Geogrid_File.py -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "----\n",
    "#### Execute the script using command-line syntax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Command to run:\n",
      "\n",
      "python Build_GeoTiff_From_Geogrid_File.py \\\n",
      "\t -i /home/docker/GIS_Training/Croton_Lambert/geo_em.d01.nc \\\n",
      "\t -v HGT_M \\\n",
      "\t -o /home/docker/GIS_Training/Outputs/HGT_M.tif\n",
      "\n",
      "Script initiated at Tue Nov  3 00:15:08 2020\n",
      "Using default variable name: HGT_M\n",
      "Input WPS Geogrid or Fulldom file: /home/docker/GIS_Training/Croton_Lambert/geo_em.d01.nc\n",
      "Input netCDF variable name: HGT_M\n",
      "Output raster file: /home/docker/GIS_Training/Outputs/HGT_M.tif\n",
      "WPS netCDF projection identification initiated...\n",
      "    Map Projection: Lambert Conformal Conic\n",
      "    Using MOAD_CEN_LAT for latitude of origin.\n",
      "    Using Standard Parallel 2 in Lambert Conformal Conic map projection.\n",
      "    Geo-referencing step completed without error in  0.05 seconds.\n",
      "    X-dimension: 'west_east'.\n",
      "    Y-dimension: 'south_north'.\n",
      "    Reversing order of dimension 'south_north'\n",
      "    Time dimension found: 'Time'.\n",
      "      Time dimension size = 1.\n",
      "    Dimensions and indices or slices on those dimensions:\n",
      "        Time: 0\n",
      "        south_north: slice(None, None, -1)\n",
      "        west_east: slice(None, None, None)\n",
      "    Size of array being sent to raster: (16, 15)\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "    Bands in output raster: 1\n",
      "    Created GTiff format raster from HGT_M variable: /home/docker/GIS_Training/Outputs/HGT_M.tif\n",
      "Process complted in 0.08 seconds.\n"
     ]
    }
   ],
   "source": [
    "# Define the variable to export to raster\n",
    "in_var = \"HGT_M\"\n",
    "\n",
    "# Define the output raster file using variable name defined above\n",
    "out_file = os.path.join(output_folder, f'{in_var}.tif')\n",
    "\n",
    "# Print information to screen for reference\n",
    "print('Command to run:\\n')\n",
    "print('python Build_GeoTiff_From_Geogrid_File.py \\\\\\n\\t -i {0} \\\\\\n\\t -v {1} \\\\\\n\\t -o {2}\\n'.format(in_geogrid, in_var, out_file))\n",
    "\n",
    "# Run the script with required parameters\n",
    "! python Build_GeoTiff_From_Geogrid_File.py -i {in_geogrid} -v {in_var} -o {out_file}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Above, you will see the elevation grid applied to the map. This grid is 1km, so there is not much detail, but it is still useful to see if the topographic features are in the correct geographic locations according to the basemap. Also, there is no color-ramp for reference. This is a limitation of using the web browser over a GIS application."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Building the hydrologic routing grids, aka the \"routing stack\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The `Build_Routing_Stack.py` script is a program to build the full set of hydrologically-processed routing grids and additional data required by WRF-Hydro. This is the main utility for performing WRF-Hydro GIS pre-processing. The required inputs are the domain file (WPS geogrid file), desired routing grid resolution as a function of geogrid resolution, and other options and parameter values. The output will be a \"routing stack\" zip file with WRF-Hydro domain and parameter files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "• Required Parameters:<br>\n",
    "&emsp;`-i` -WRF/WPS GEOGRID file (geo_em.d0*.nc)<br>\n",
    "&emsp;`-d` -High-resolution Elevation raster file (Esri GRID, GeoTIFF, VRT, etc.)<br>\n",
    "&emsp;`-R` -Regridding Factor – nesting relationship of routing:land surface model grid cells<br>\n",
    "&emsp;`-t` -Minimum basin area threshold (in routing grid cells)<br>\n",
    "&emsp;`-o` -Output ZIP File containing all script outputs<br>\n",
    "\n",
    "• Optional Parameters:<br>\n",
    "&emsp;`--CSV` -Station Locations location file (.csv)<br>\n",
    "&emsp;`-b` -Option to mask channel grids not contributing to provided station locations<br>\n",
    "&emsp;`-r` -Reach based (Muskingum / Muskingum-Cunge) routing option<br>\n",
    "&emsp;`-l` -Lake Polygons (polygon feature class or .shp)<br>\n",
    "&emsp;`-O` -OVROUGHRTFAC – Multiplier on Manning's roughness for overland flow. default=1.0<br>\n",
    "&emsp;`-T` -RETDEPRTFAC – Multiplier on maximum retention depth before flow is routed as overland flow. default=1.0<br>\n",
    "&emsp;-LKSATFAC – (script global variable) Multiplier on saturated hydraulic conductivity in lateral flow direction. default=1000.0<br>\n",
    "&emsp;`--starts` -Path to point shapefile or feature class containing channel initiation locations (overrides `-t` parameter)<br>\n",
    "&emsp;`--gw` -Path to polygon shapefile or feature class containing prescribed groundwater basin locations<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Request help message from the script\n",
    "This tool has many different parameters and possible configurations. Using the command line, we can take a look at the different arguements that can be used, if they are required or optional, and what their default values are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Script initiated at Tue Nov  3 00:15:09 2020\n",
      "usage: Build_Routing_Stack.py [-h] -i IN_GEOGRID [--CSV IN_CSV]\n",
      "                              [-b BASIN_MASK] [-r RB_ROUTING]\n",
      "                              [-l IN_RESERVOIRS] -d INDEM [-R CELLSIZE]\n",
      "                              [-t THRESHOLD] [-o OUT_ZIP_FILE]\n",
      "                              [-O OVROUGHRTFAC_VAL] [-T RETDEPRTFAC_VAL]\n",
      "                              [--starts CHANNEL_STARTS] [--gw GW_POLYS]\n",
      "\n",
      "This is a program to perform the full routing-stack GIS pre-processingfor WRF-\n",
      "Hydro. The inputs will be related to the domain, the desired routing nest\n",
      "factor, and other options and parameter values. The output will be a routing\n",
      "stack zip file with WRF-Hydro domain and parameter files.\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help            show this help message and exit\n",
      "  -i IN_GEOGRID         Path to WPS geogrid (geo_em.d0*.nc) file [REQUIRED]\n",
      "  --CSV IN_CSV          Path to input forecast point CSV file [OPTIONAL]\n",
      "  -b BASIN_MASK         Mask CHANNELGRID variable to forecast basins?\n",
      "                        [True/False]. default=False\n",
      "  -r RB_ROUTING         Create reach-based routing (RouteLink) files?\n",
      "                        [True/False]. default=False\n",
      "  -l IN_RESERVOIRS      Path to reservoirs shapefile or feature class\n",
      "                        [OPTIONAL]. If -l is TRUE, this is required.\n",
      "  -d INDEM              Path to input high-resolution elevation raster\n",
      "                        [REQUIRED]\n",
      "  -R CELLSIZE           Regridding (nest) Factor. default=10\n",
      "  -t THRESHOLD          Number of routing grid cells to define stream.\n",
      "                        default=200\n",
      "  -o OUT_ZIP_FILE       Output routing stack ZIP file\n",
      "  -O OVROUGHRTFAC_VAL   OVROUGHRTFAC value. default=1.0\n",
      "  -T RETDEPRTFAC_VAL    RETDEPRTFAC value. default=1.0\n",
      "  --starts CHANNEL_STARTS\n",
      "                        Path to channel initiation points feature class. Must\n",
      "                        be 2D point type. [OPTIONAL]\n",
      "  --gw GW_POLYS         Path to groundwater polygons feature class [OPTIONAL]\n"
     ]
    }
   ],
   "source": [
    "! python Build_Routing_Stack.py -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "#### Execute the script using command-line syntax\n",
    "We will begin by assigning file paths to python variables, then substitude those values in the command line syntax."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Command to run:\n",
      "\n",
      "python Build_Routing_Stack.py \\\n",
      "\t -i /home/docker/GIS_Training/Croton_Lambert/geo_em.d01.nc \\\n",
      "\t -l /home/docker/GIS_Training/Croton_Lambert/lake_shapes/lakes.shp \\\n",
      "\t --CSV /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv \\\n",
      "\t -d /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif \\\n",
      "\t -R 4 \\\n",
      "\t -t 20 \\\n",
      "\t -o /home/docker/GIS_Training/Outputs/croton_test.zip\n",
      "\n",
      "Script initiated at Tue Nov  3 00:15:11 2020\n",
      "  Parameter values that have not been altered from script default values:\n",
      "    Using default basin mask setting: False\n",
      "    Using default reach-based routing setting: False\n",
      "    Using default OVROUGHRTFAC parameter value: 1.0\n",
      "    Using default RETDEPRTFAC parameter value: 1.0\n",
      "  Values that will be used in building this routing stack:\n",
      "    Input WPS Geogrid file: /home/docker/GIS_Training/Croton_Lambert/geo_em.d01.nc\n",
      "    Forecast Point CSV file: /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv\n",
      "    Mask CHANNELGRID variable to forecast basins?: False\n",
      "    Create reach-based routing (RouteLink) files?: False\n",
      "    Lake polygon feature class: /home/docker/GIS_Training/Croton_Lambert/lake_shapes/lakes.shp\n",
      "    Input high-resolution DEM: /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif\n",
      "    Regridding factor: 4\n",
      "    Stream initiation threshold: 20\n",
      "    OVROUGHRTFAC parameter value: 1.0\n",
      "    RETDEPRTFAC parameter value: 1.0\n",
      "    Input channel initiation start point feature class: None\n",
      "    Input groundwater basin polygons: None\n",
      "    Output ZIP file: /home/docker/GIS_Training/Outputs/croton_test.zip\n",
      "  Running Process GEOGRID function\n",
      "    Forecast points provided.\n",
      "  Reach-based routing files will not be created.\n",
      "WPS netCDF projection identification initiated...\n",
      "    Map Projection: Lambert Conformal Conic\n",
      "    Using MOAD_CEN_LAT for latitude of origin.\n",
      "    Using Standard Parallel 2 in Lambert Conformal Conic map projection.\n",
      "    Geo-referencing step completed without error in  0.04 seconds.\n",
      "  Building sub-grid of model grid.\n",
      "    Original grid spacing dx=1000.0, dy=-1000.0\n",
      "    Original grid size: rows=16, cols=15\n",
      "    New grid spacing: dx=250.0, dy=-250.0\n",
      "    New dimensions: rows=64, cols=60\n",
      "    Created projection definition from input NetCDF GEOGRID file.\n",
      "    Proj4: +proj=lcc +lat_0=41.4710083007812 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs\n",
      "    Coarse grid GeoTransform: 1841999.0194101864 1000.0 0 278495.814356732 0 -1000.0\n",
      "    Coarse grid extent [Xmin, Ymin, Xmax, Ymax]: [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732]\n",
      "    Fine grid extent [Xmin, Ymin, Xmax, Ymax]:   [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732]\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "  Creating CF-netCDF File.\n",
      "    Esri PE String: PROJCS[\"unnamed\",GEOGCS[\"GCS_Sphere\",DATUM[\"D_Sphere\",SPHEROID[\"unnamed\",6370000.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-97.0],PARAMETER[\"Standard_Parallel_1\",30.0],PARAMETER[\"Standard_Parallel_2\",60.0],PARAMETER[\"Latitude_Of_Origin\",41.4710083007812],UNIT[\"Meter\",1.0]]\n",
      "    Map Projection of input raster : lambert_conformal_conic\n",
      "    Starting Process: Building to XMap/YMap\n",
      "    Conversion of input raster to XMap/YMap completed without error.\n",
      "  netCDF global attributes set after  0.01 seconds.\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.13 seconds.\n",
      "  Deriving geocentric coordinates on routing grid from bilinear interpolation of geogrid coordinates.\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "  Creating CF-netCDF File.\n",
      "    Esri PE String: PROJCS[\"unnamed\",GEOGCS[\"GCS_Sphere\",DATUM[\"D_Sphere\",SPHEROID[\"unnamed\",6370000.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-97.0],PARAMETER[\"Standard_Parallel_1\",30.0],PARAMETER[\"Standard_Parallel_2\",60.0],PARAMETER[\"Latitude_Of_Origin\",41.4710083007812],UNIT[\"Meter\",1.0]]\n",
      "    Map Projection of input raster : lambert_conformal_conic\n",
      "    Starting Process: Building to XMap/YMap\n",
      "    Conversion of input raster to XMap/YMap completed without error.\n",
      "    Proceeding to add LATITUDE and LONGITUDE variables after     0.01 seconds.\n",
      "  netCDF global attributes set after  0.01 seconds.\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "    Process: landuse written to output netCDF.\n",
      "Terrain processing step initiated...\n",
      "    Using WhiteboxTools v1.4.0 by Dr. John B. Lindsay (c) 2017-2020\n",
      "    Depression Filling algorithm: Whitebox Fill Depressions.\n",
      "    Process: TOPOGRAPHY written to output netCDF.\n",
      "    Coerced 21 0-value flow direction cells to flow off of the grid.\n",
      "    Process: FLOWDIRECTION written to output netCDF.\n",
      "    Process: FLOWACC written to output netCDF.\n",
      "    Flow accumulation will be thresholded to build channel pixels.\n",
      "    Process: CHANNELGRID written to output netCDF.\n",
      "    Process: STREAMORDER written to output netCDF.\n",
      "    Process: RETDEPRTFAC written to output netCDF.\n",
      "    Process: OVROUGHRTFAC written to output netCDF.\n",
      "    Process: LKSATFAC written to output netCDF.\n",
      "Terrain processing step completed without error in  0.05 seconds.\n",
      "    Forecast points provided and basins being delineated.\n",
      "    Process: frxst_pts written to output netCDF.\n",
      "    Process: basn_msk written to output netCDF.\n",
      "    Channelgrid will not be masked to basins.\n",
      "    Built forecast point outputs in  0.11 seconds.\n",
      "    Reservoir polygons provided. Lake routing will be activated.\n",
      "      Adding reservoirs to routing stack.\n",
      "      Gridded: True\n",
      "      MEMORY driver is available.\n",
      "  Done producing output vector polygon shapefile in  0.00 seconds\n",
      "    Input shapefile projection does not match requested output. Transforming.\n",
      "    Number of output features: 1 of 1\n",
      "  Completed reprojection and-or clipping in 0.16 seconds.\n",
      "    Adding auto-incremented lake ID field (1...n)\n",
      "    Starting to gather lake centroid and area information.\n",
      "    Done gathering lake centroid information.\n",
      "    Found 1 lakes on active channels. Lost 0 lakes that were not on active channels.\n",
      "    Removing lakes not on gridded channel network\n",
      "    Process: LAKEGRID written to output netCDF.\n",
      "    Process: CHANNELGRID written to output netCDF.\n",
      "    Starting to create lake parameter table.\n",
      "        Lakes Table: 1 Lakes\n",
      "        Starting to fill in lake parameter table NC file.\n",
      "        Done writing LAKEPARM.nc table to disk.\n",
      "    Lake parameter table created without error in  0.31 seconds.\n",
      "Beginning to build 2D groundwater basin inputs\n",
      "  Building groundwater inputs using FullDom LINKID local basins\n",
      "    Generating LINKID grid for building local sub-basins.\n",
      "Finished building fine-grid groundwater basin grids in  0.02 seconds\n",
      "Beginning to build coarse-grid groundwater basins and parameters\n",
      "    Found 72 basins in the watershed grid\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 1000.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "    Found 69 basins (potentially including nodata values) in the file after resampling to the coarse grid.\n",
      "    GDAL Data type derived from input array: 5 (int32)\n",
      "  Creating CF-netCDF File.\n",
      "    Esri PE String: PROJCS[\"unnamed\",GEOGCS[\"GCS_Sphere\",DATUM[\"D_Sphere\",SPHEROID[\"unnamed\",6370000.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-97.0],PARAMETER[\"Standard_Parallel_1\",30.0],PARAMETER[\"Standard_Parallel_2\",60.0],PARAMETER[\"Latitude_Of_Origin\",41.4710083007812],UNIT[\"Meter\",1.0]]\n",
      "    Map Projection of input raster : lambert_conformal_conic\n",
      "    Starting Process: Building to XMap/YMap\n",
      "    Conversion of input raster to XMap/YMap completed without error.\n",
      "  netCDF global attributes set after  0.01 seconds.\n",
      "    NC dimensions: 16, 15\n",
      "    GWBUCKS array dimensions: 16, 15\n",
      "    Process: /home/docker/GIS_Training/Outputs/scratchdir/GWBASINS.nc completed without error\n",
      "    Finished building groundwater grid file in  0.01 seconds\n",
      "    Calculating size and ID parameters for basin polygons.\n",
      "    Created output bucket parameter table (.nc): /home/docker/GIS_Training/Outputs/scratchdir/GWBUCKPARM.nc.\n",
      "  Finished building groundwater bucket parameter table in  0.00 seconds.\n",
      "Finished building groundwater parameter files in  0.03 seconds\n",
      "Built output .zip file in  0.81 seconds.\n",
      "Process completed in 0.81 seconds.\n"
     ]
    }
   ],
   "source": [
    "import Build_Routing_Stack\n",
    "\n",
    "# Define script input parameters using python variables\n",
    "in_geogrid = os.path.join(data_folder, 'geo_em.d01.nc')\n",
    "lakes = os.path.join(data_folder, 'lake_shapes', 'lakes.shp')\n",
    "csv = os.path.join(data_folder, 'croton_frxst_pts_FOSS.csv')\n",
    "in_dem = os.path.join(data_folder, 'NED_30m_Croton.tif')\n",
    "regrid_factor = 4\n",
    "routing_cells = 20\n",
    "out_zip = os.path.join(output_folder, 'croton_test.zip')\n",
    "\n",
    "# Print information to screen for reference\n",
    "print('Command to run:\\n')\n",
    "print('python Build_Routing_Stack.py \\\\\\n\\t -i {0} \\\\\\n\\t -l {1} \\\\\\n\\t --CSV {2} \\\\\\n\\t -d {3} \\\\\\n\\t -R {4} \\\\\\n\\t -t {5} \\\\\\n\\t -o {6}\\n'.format(in_geogrid, lakes, csv, in_dem, regrid_factor, routing_cells, out_zip))\n",
    "\n",
    "# Run the script with required parameters\n",
    "! python Build_Routing_Stack.py -i {in_geogrid} -l {lakes} --CSV {csv} -d {in_dem} -R {regrid_factor} -t {routing_cells} -o {out_zip}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Understanding the outputs\n",
    "\n",
    "The `Build_Routing_Stack.py` script creates a Zip archive of output files according to the options provided to the tool. There will be at least four netCDF files. The output Zip file may additionally include shapefiles (.shp and accompanying files) describing the geometry of modeled lakes or the vector stream network. Below is an alphabetically sorted list of gridded variables that are created by the `Build_Routing_Stack.py` tool.\n",
    "\n",
    "Fulldom_hires.nc file. This file stores 2D gridded variables that describe the hydro routing grid:\n",
    "  + CHANNELGRID - The channel grid. Channel pixels = 0, non-channel pixels = -9999. If the `-b` option is set to TRUE, the output will be masked to the gaged basins provided, where non-gaged channels are given a value of ‘-1’. If lake routing is activated, lake outflow points will be identified by the lake ID value.\n",
    "  + FLOWACC – Flow accumulation grid. This grid gives the number of contributing cells for each cell in the domain. This grid is provided for convenience and is not read by WRF-Hydro.\n",
    "  + FLOWDIRECTION – Flow direction grid. This grid gives the direction of flow using the D8 algorithm between each cell and the steepest downslope neighbor according to Jenson and Domingue (1988). The result is an integer grid with values ranging from 1 to 128.\n",
    "  + frxst_pts – Gage location grid. If a forecast point CSV file is provided, the grid will have a cell identified at the location of each forecast point (gage) in the gage CSV file. If no input CSV gage location file is provided, this grid will be uniform with values of ‘-9999’. Gage pixels are numbered in the same way as the ‘basn_msk’ grid. NoData cells are given a value of ‘-9999’.\n",
    "  + basn_msk – Forecast basins grid. If a CSV gage location file is provided, catchments are delineated from a point that is up to 3 pixels downstream of the gage coordinates. This distance can be modified by altering the ‘walker’ global variable in the ‘wrfhydro_functions.py’ script. If masking of the ‘CHANNELGRID’ is selected, this layer is the mask. Basins are numbered according to the values in the ‘FID’ field of the input gage CSV file. If no gage location file is provided, this grid will be uniform with values of ‘-9999’.\n",
    "  + LAKEGRID – The lake grid. If a lake polygon shapefile is provided to the `-l` parameter, this grid will contain ID values for each lake that can be resolved on the routing grid. Otherwise, this grid will be uniform with values of -9999.\n",
    "  + landuse – This is the same data as the ‘LU_INDEX’ variable in the GEOGRID file, but resampled using Nearest Neighbor assignment to the resolution of the routing grid. This grid is provided for convenience and is not read by WRF-Hydro.\n",
    "  + LATITUDE – Grid of the latitude at the center of each grid cell, in a geographic coordinate system (WGS84).\n",
    "  + LINKID – The channel ID grid. This grid provides a unique integer identifier for each channel segment that is defined in the ‘link’ variable of the ‘Route_Link.nc’ file and the ‘STRM_VAL’ field in the ‘streams.shp’ shapefile. The ‘LINKID’ grid will only be created if the option `-r` is TRUE.\n",
    "  + LONGITUDE – Grid of longitude value at the center of each grid cell, in a geographic coordinate system (WGS84).\n",
    "  + OVROUGHRTFAC – OVROUGHRTFAC parameter. Currently set to a default of 1.0. This default value may be changed by providing an alternate value to the `-O` parameter.\n",
    "  + RETDEPRTFAC – RETDEPRTFAC (retention depth multiplier) parameter. Currently set to a default of 1.0. This default value may be changed by providing an alternate value `-T` parameter.\n",
    "  + STREAMORDER – Stream order grid, calculated using the Strahler method (Strahler 1957).\n",
    "  + TOPOGRAPHY – Elevation grid. The units of elevation are the same as the input elevation raster dataset (`-d INDEM`), which should be specified in meters (m) above sea level (ASL). This grid is derived from the elevation values in the input elevation raster, but has been resampled to the routing grid resolution, and pit-filled to remove depressions.\n",
    "\n",
    "Other files:\n",
    "  + GEOGRID_LDASOUT_Spatial_Metadata.nc  - This is a CF-netCDF format file that provides the spatial metadata associated with the GEOGRID variables, which contsitute the LSM grid. By default, no 2-dimensional grids are written to the file. This file may be used by WRF-Hydro for appending geospatial metadata to the land surface model output, if necessary.\n",
    "  + GWBASINS.nc - This is a 2D netCDF file of the location of groundwater basins, regridded to the LSM grid resolution. NoData cells are given a value of ‘-9999’. This file is by default created using the ‘FullDom LINKID local basins’ method of defining groundwater basins.\n",
    "  + GWBUCKPARM.nc - The 1D groundwater basin parameter file.\n",
    "  + LAKEPARM.nc – Lake parameter table. This 1D netCDF format file is created if a lake shapefile is provided as input to the `-l` parameter. The table will contain a record for each lake feature in the + Fulldom_hires.nc ‘LAKEGRID’ variable, and contain derived and default parameters for each lake.\n",
    "  + Route_Link.nc – The reach-based routing parameter file. This 1D netCDF format file is created if the `-r` parameter is TRUE. The file contains a record for each stream segment. The stream segments in this table are also identified by the unique ‘LINKID’ values in the ‘LINKID’ variable in the ‘Fulldom_hires.nc’ file, and values in the ‘STRM_VAL’ field of the output ‘streams.shp’ shapefile. This table contains derived and default stream segment parameters that are calculated based on the vector stream network and topology in the ‘streams.shp’ shapefile.\n",
    "  + streams.* (ancillary) - Streams shapefile, containing one feature for each stream segment in the domain. This file is meant to accompany the ‘Route_Link.nc’ reach-based routing parameter file and Fulldom_hires.nc ‘LINKID’ variable. The ‘streams’ shapefile is only created when the option `-l` is used. The ‘STRM_VAL’ field is the unique identifier for each stream segment, and corresponds to the ‘link’ variable of the ‘Route_Link.nc’ file and the ‘LINKID’ variable in the ‘Fulldom_hires.nc’ file. The geometry of the stream segments in this shapefile informs many of the parameters in the ‘Route_Link.nc’ file.\n",
    "  + lakes.* (ancillary) - Lakes shapefile, containing one feature for each reservoir in the simulation domain. This file is meant to accopmany the 'LAKEPARM.nc' reservoir parameter file. If `-r TRUE` is used, then Fulldom_hires.nc 'LAKEID' variable will contain -9999 values only. The geometry of reservoir objects informs many of the parameters in 'LAKEPARM.nc' file.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Examine outputs of GIS pre-processor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The `Examine_Outputs_of_GIS_Preprocessor.py` script is a tool used to create easy-to-use geospatial files for examining the resulting WRF-Hydro input files in a GIS. This script takes the output ZIP file generated using the Process Geogrid script (executed above) and creates a raster from each 2D variable in each WRF-Hydro input netCDF file. In addition, other data will be extracted such as any shapefiles, 1D netCDF tables, etc. The input to the tool should be a .zip file that was created using the WRF Hydro pre-processing tools. The tool will create the output folder if it does not already exist, and write all results to that location.\n",
    "\n",
    "#### Request help message from the script\n",
    "This tool has a single input and single output parameter. Using the command line, we can take a look at the arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Script initiated at Tue Nov  3 00:15:13 2020\n",
      "usage: Examine_Outputs_of_GIS_Preprocessor.py [-h] -i IN_ZIP -o OUT_FOLDER\n",
      "\n",
      "This tool takes the output zip file from the ProcessGeogrid script and creates\n",
      "a raster from each output NetCDF file. The Input should be a .zip file that\n",
      "was created using theWRF Hydro pre-processing tools. The tool will create the\n",
      "folder which will contain theresults (out_folder), if that folder does not\n",
      "already exist.\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help     show this help message and exit\n",
      "  -i IN_ZIP      Path to WRF Hydro routing grids zip file.\n",
      "  -o OUT_FOLDER  Path to output folder.\n"
     ]
    }
   ],
   "source": [
    "! python Examine_Outputs_of_GIS_Preprocessor.py -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "#### Execute the script using command-line syntax\n",
    "We will define the output directory for the tool to create, then use defined variables to execute the command line tool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Command to run:\n",
      "\n",
      "python Examine_Outputs_of_GIS_Preprocessor.py \\\n",
      "\t -i /home/docker/GIS_Training/Outputs/croton_test.zip \\\n",
      "\t -o /home/docker/GIS_Training/Outputs/Raster_Outputs\n",
      "\n",
      "Script initiated at Tue Nov  3 00:15:14 2020\n",
      "    File Copied: lakes.prj\n",
      "    File Copied: lakes.shx\n",
      "  GeoTransform: (1841999.0194101864, 1000.0, 0.0, 278495.814356732, 0.0, -1000.0)\n",
      "  DX: 1000.0\n",
      "  DY: 1000.0\n",
      "  PROJ.4 string: +proj=lcc +lat_0=41.4710083007812 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/BASIN.tif\n",
      "    File Copied: LAKEPARM.nc\n",
      "    File Copied: lakes.dbf\n",
      "    File Copied: GEOGRID_LDASOUT_Spatial_Metadata.nc\n",
      "    File Copied: lakes.shp\n",
      "  GeoTransform: (1841999.0194101864, 250.0, 0.0, 278495.814356732, 0.0, -250.0)\n",
      "  DX: 250.0\n",
      "  DY: 250.0\n",
      "  PROJ.4 string: +proj=lcc +lat_0=41.4710083007812 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/CHANNELGRID.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/FLOWDIRECTION.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/FLOWACC.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/TOPOGRAPHY.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/RETDEPRTFAC.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/OVROUGHRTFAC.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/STREAMORDER.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/frxst_pts.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/basn_msk.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/LAKEGRID.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/landuse.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/LKSATFAC.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/LATITUDE.tif\n",
      "    File Created: /home/docker/GIS_Training/Outputs/Raster_Outputs/LONGITUDE.tif\n",
      "    File Copied: GWBUCKPARM.nc\n",
      "Extraction of WRF routing grids completed.\n",
      "Process complted in 0.10 seconds.\n"
     ]
    }
   ],
   "source": [
    "# Define output directory to store GeoTiff output of all routing stack grids\n",
    "raster_outputs = os.path.join(output_folder, \"Raster_Outputs\")\n",
    "\n",
    "# Print information to screen for reference\n",
    "print('Command to run:\\n')\n",
    "print('python Examine_Outputs_of_GIS_Preprocessor.py \\\\\\n\\t -i {0} \\\\\\n\\t -o {1}\\n'.format(out_zip, raster_outputs))\n",
    "\n",
    "# Run the script with required parameters\n",
    "! python Examine_Outputs_of_GIS_Preprocessor.py -i {out_zip} -o {raster_outputs}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. [Optional] Build Non-NWM WRF-Hydro Configurations of the Croton, NY Test Case\n",
    "\n",
    "We will run command-line arguments to create output directories for each of the test cases. Then, each cell will run the WRF-Hydro GIS Pre-processor to build the 'routing stack' files. The final command is used to unzip the outputs so they can be more readily used by WRF-Hydro. There are 4 Non-NWM test-case configurations which can all be created using different combinations of arguments to the GIS pre-processor. Currently missing at this time is the Reach-with-Lakes capability as well as the ability to build wrfinput, hydro2dtbl, and soil_properties files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A. Gridded channel routing configuration with reservoirs and forecast points\n",
    "\n",
    "This domain is a gridded domain with lakes, gages and a regridding factor of 4 and threshold of 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Script initiated at Tue Nov  3 00:15:17 2020\n",
      "  Parameter values that have not been altered from script default values:\n",
      "    Using default basin mask setting: False\n",
      "    Using default reach-based routing setting: False\n",
      "    Using default OVROUGHRTFAC parameter value: 1.0\n",
      "    Using default RETDEPRTFAC parameter value: 1.0\n",
      "  Values that will be used in building this routing stack:\n",
      "    Input WPS Geogrid file: /home/docker/wrf-hydro-training/example_case/Gridded/DOMAIN/geo_em.d01.nc\n",
      "    Forecast Point CSV file: /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv\n",
      "    Mask CHANNELGRID variable to forecast basins?: False\n",
      "    Create reach-based routing (RouteLink) files?: False\n",
      "    Lake polygon feature class: /home/docker/GIS_Training/Croton_Lambert/lake_shapes/lakes.shp\n",
      "    Input high-resolution DEM: /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif\n",
      "    Regridding factor: 4\n",
      "    Stream initiation threshold: 20\n",
      "    OVROUGHRTFAC parameter value: 1.0\n",
      "    RETDEPRTFAC parameter value: 1.0\n",
      "    Input channel initiation start point feature class: None\n",
      "    Input groundwater basin polygons: None\n",
      "    Output ZIP file: /home/docker/GIS_Training/Outputs/Gridded/Gridded_r4_t20_lakes_frxst_mask.zip\n",
      "  Running Process GEOGRID function\n",
      "    Forecast points provided.\n",
      "  Reach-based routing files will not be created.\n",
      "WPS netCDF projection identification initiated...\n",
      "    Map Projection: Lambert Conformal Conic\n",
      "    Using MOAD_CEN_LAT for latitude of origin.\n",
      "    Using Standard Parallel 2 in Lambert Conformal Conic map projection.\n",
      "    Geo-referencing step completed without error in  0.05 seconds.\n",
      "  Building sub-grid of model grid.\n",
      "    Original grid spacing dx=1000.0, dy=-1000.0\n",
      "    Original grid size: rows=16, cols=15\n",
      "    New grid spacing: dx=250.0, dy=-250.0\n",
      "    New dimensions: rows=64, cols=60\n",
      "    Created projection definition from input NetCDF GEOGRID file.\n",
      "    Proj4: +proj=lcc +lat_0=41.4710083007812 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs\n",
      "    Coarse grid GeoTransform: 1841999.0194101864 1000.0 0 278495.814356732 0 -1000.0\n",
      "    Coarse grid extent [Xmin, Ymin, Xmax, Ymax]: [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732]\n",
      "    Fine grid extent [Xmin, Ymin, Xmax, Ymax]:   [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732]\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "  Creating CF-netCDF File.\n",
      "    Esri PE String: PROJCS[\"unnamed\",GEOGCS[\"GCS_Sphere\",DATUM[\"D_Sphere\",SPHEROID[\"unnamed\",6370000.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-97.0],PARAMETER[\"Standard_Parallel_1\",30.0],PARAMETER[\"Standard_Parallel_2\",60.0],PARAMETER[\"Latitude_Of_Origin\",41.4710083007812],UNIT[\"Meter\",1.0]]\n",
      "    Map Projection of input raster : lambert_conformal_conic\n",
      "    Starting Process: Building to XMap/YMap\n",
      "    Conversion of input raster to XMap/YMap completed without error.\n",
      "  netCDF global attributes set after  0.01 seconds.\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.13 seconds.\n",
      "  Deriving geocentric coordinates on routing grid from bilinear interpolation of geogrid coordinates.\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "  Creating CF-netCDF File.\n",
      "    Esri PE String: PROJCS[\"unnamed\",GEOGCS[\"GCS_Sphere\",DATUM[\"D_Sphere\",SPHEROID[\"unnamed\",6370000.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-97.0],PARAMETER[\"Standard_Parallel_1\",30.0],PARAMETER[\"Standard_Parallel_2\",60.0],PARAMETER[\"Latitude_Of_Origin\",41.4710083007812],UNIT[\"Meter\",1.0]]\n",
      "    Map Projection of input raster : lambert_conformal_conic\n",
      "    Starting Process: Building to XMap/YMap\n",
      "    Conversion of input raster to XMap/YMap completed without error.\n",
      "    Proceeding to add LATITUDE and LONGITUDE variables after     0.01 seconds.\n",
      "  netCDF global attributes set after  0.01 seconds.\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "    Process: landuse written to output netCDF.\n",
      "Terrain processing step initiated...\n",
      "    Using WhiteboxTools v1.4.0 by Dr. John B. Lindsay (c) 2017-2020\n",
      "    Depression Filling algorithm: Whitebox Fill Depressions.\n",
      "    Process: TOPOGRAPHY written to output netCDF.\n",
      "    Coerced 21 0-value flow direction cells to flow off of the grid.\n",
      "    Process: FLOWDIRECTION written to output netCDF.\n",
      "    Process: FLOWACC written to output netCDF.\n",
      "    Flow accumulation will be thresholded to build channel pixels.\n",
      "    Process: CHANNELGRID written to output netCDF.\n",
      "    Process: STREAMORDER written to output netCDF.\n",
      "    Process: RETDEPRTFAC written to output netCDF.\n",
      "    Process: OVROUGHRTFAC written to output netCDF.\n",
      "    Process: LKSATFAC written to output netCDF.\n",
      "Terrain processing step completed without error in  0.04 seconds.\n",
      "    Forecast points provided and basins being delineated.\n",
      "    Process: frxst_pts written to output netCDF.\n",
      "    Process: basn_msk written to output netCDF.\n",
      "    Channelgrid will not be masked to basins.\n",
      "    Built forecast point outputs in  0.11 seconds.\n",
      "    Reservoir polygons provided. Lake routing will be activated.\n",
      "      Adding reservoirs to routing stack.\n",
      "      Gridded: True\n",
      "      MEMORY driver is available.\n",
      "  Done producing output vector polygon shapefile in  0.00 seconds\n",
      "    Input shapefile projection does not match requested output. Transforming.\n",
      "    Number of output features: 1 of 1\n",
      "  Completed reprojection and-or clipping in 0.16 seconds.\n",
      "    Adding auto-incremented lake ID field (1...n)\n",
      "    Starting to gather lake centroid and area information.\n",
      "    Done gathering lake centroid information.\n",
      "    Found 1 lakes on active channels. Lost 0 lakes that were not on active channels.\n",
      "    Removing lakes not on gridded channel network\n",
      "    Process: LAKEGRID written to output netCDF.\n",
      "    Process: CHANNELGRID written to output netCDF.\n",
      "    Starting to create lake parameter table.\n",
      "        Lakes Table: 1 Lakes\n",
      "        Starting to fill in lake parameter table NC file.\n",
      "        Done writing LAKEPARM.nc table to disk.\n",
      "    Lake parameter table created without error in  0.29 seconds.\n",
      "Beginning to build 2D groundwater basin inputs\n",
      "  Building groundwater inputs using FullDom LINKID local basins\n",
      "    Generating LINKID grid for building local sub-basins.\n",
      "Finished building fine-grid groundwater basin grids in  0.02 seconds\n",
      "Beginning to build coarse-grid groundwater basins and parameters\n",
      "    Found 72 basins in the watershed grid\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 1000.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "    Found 69 basins (potentially including nodata values) in the file after resampling to the coarse grid.\n",
      "    GDAL Data type derived from input array: 5 (int32)\n",
      "  Creating CF-netCDF File.\n",
      "    Esri PE String: PROJCS[\"unnamed\",GEOGCS[\"GCS_Sphere\",DATUM[\"D_Sphere\",SPHEROID[\"unnamed\",6370000.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-97.0],PARAMETER[\"Standard_Parallel_1\",30.0],PARAMETER[\"Standard_Parallel_2\",60.0],PARAMETER[\"Latitude_Of_Origin\",41.4710083007812],UNIT[\"Meter\",1.0]]\n",
      "    Map Projection of input raster : lambert_conformal_conic\n",
      "    Starting Process: Building to XMap/YMap\n",
      "    Conversion of input raster to XMap/YMap completed without error.\n",
      "  netCDF global attributes set after  0.00 seconds.\n",
      "    NC dimensions: 16, 15\n",
      "    GWBUCKS array dimensions: 16, 15\n",
      "    Process: /home/docker/GIS_Training/Outputs/Gridded/scratchdir/GWBASINS.nc completed without error\n",
      "    Finished building groundwater grid file in  0.01 seconds\n",
      "    Calculating size and ID parameters for basin polygons.\n",
      "    Created output bucket parameter table (.nc): /home/docker/GIS_Training/Outputs/Gridded/scratchdir/GWBUCKPARM.nc.\n",
      "  Finished building groundwater bucket parameter table in  0.00 seconds.\n",
      "Finished building groundwater parameter files in  0.03 seconds\n",
      "Built output .zip file in  0.77 seconds.\n",
      "Process completed in 0.77 seconds.\n",
      "Archive:  /home/docker/GIS_Training/Outputs/Gridded/Gridded_r4_t20_lakes_frxst_mask.zip\n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded/lakes.prj  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded/lakes.shx  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded/GWBASINS.nc  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded/LAKEPARM.nc  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded/lakes.dbf  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded/GEOGRID_LDASOUT_Spatial_Metadata.nc  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded/lakes.shp  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded/Fulldom_hires.nc  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded/GWBUCKPARM.nc  \n"
     ]
    }
   ],
   "source": [
    "# Create the output directory and ensure it is empty\n",
    "! mkdir -p /home/docker/GIS_Training/Outputs/Gridded\n",
    "! rm -rf /home/docker/GIS_Training/Outputs/Gridded/*\n",
    "\n",
    "# Run the GIS Pre-processing script (with line-breaks)\n",
    "! python Build_Routing_Stack.py \\\n",
    "    -i /home/docker/wrf-hydro-training/example_case/Gridded/DOMAIN/geo_em.d01.nc \\\n",
    "    -l /home/docker/GIS_Training/Croton_Lambert/lake_shapes/lakes.shp \\\n",
    "    --CSV /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv \\\n",
    "    -d /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif \\\n",
    "    -R 4 \\\n",
    "    -t 20 \\\n",
    "    -o /home/docker/GIS_Training/Outputs/Gridded/Gridded_r4_t20_lakes_frxst_mask.zip\n",
    "\n",
    "# Unzip the directory in-place\n",
    "! unzip /home/docker/GIS_Training/Outputs/Gridded/Gridded_r4_t20_lakes_frxst_mask.zip -d /home/docker/GIS_Training/Outputs/Gridded"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B. Gridded channel routing configuration with forecast points, and channel masking\n",
    "\n",
    "This domain is a gridded domain without lakes. It contains gages, masked basins and a regridding factor of 4 and thrshold of 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Script initiated at Tue Nov  3 00:15:41 2020\n",
      "  Parameter values that have not been altered from script default values:\n",
      "    Using default basin mask setting: False\n",
      "    Using default reach-based routing setting: False\n",
      "    Using default OVROUGHRTFAC parameter value: 1.0\n",
      "    Using default RETDEPRTFAC parameter value: 1.0\n",
      "  Values that will be used in building this routing stack:\n",
      "    Input WPS Geogrid file: /home/docker/wrf-hydro-training/example_case/Gridded_no_lakes/DOMAIN/geo_em.d01.nc\n",
      "    Forecast Point CSV file: /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv\n",
      "    Mask CHANNELGRID variable to forecast basins?: False\n",
      "    Create reach-based routing (RouteLink) files?: False\n",
      "    Lake polygon feature class: \n",
      "    Input high-resolution DEM: /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif\n",
      "    Regridding factor: 4\n",
      "    Stream initiation threshold: 20\n",
      "    OVROUGHRTFAC parameter value: 1.0\n",
      "    RETDEPRTFAC parameter value: 1.0\n",
      "    Input channel initiation start point feature class: None\n",
      "    Input groundwater basin polygons: None\n",
      "    Output ZIP file: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/Gridded_r4_t20_frxst_mask.zip\n",
      "  Running Process GEOGRID function\n",
      "    Forecast points provided.\n",
      "  Reach-based routing files will not be created.\n",
      "WPS netCDF projection identification initiated...\n",
      "    Map Projection: Lambert Conformal Conic\n",
      "    Using MOAD_CEN_LAT for latitude of origin.\n",
      "    Using Standard Parallel 2 in Lambert Conformal Conic map projection.\n",
      "    Geo-referencing step completed without error in  0.04 seconds.\n",
      "  Building sub-grid of model grid.\n",
      "    Original grid spacing dx=1000.0, dy=-1000.0\n",
      "    Original grid size: rows=16, cols=15\n",
      "    New grid spacing: dx=250.0, dy=-250.0\n",
      "    New dimensions: rows=64, cols=60\n",
      "    Created projection definition from input NetCDF GEOGRID file.\n",
      "    Proj4: +proj=lcc +lat_0=41.4710083007812 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs\n",
      "    Coarse grid GeoTransform: 1841999.0194101864 1000.0 0 278495.814356732 0 -1000.0\n",
      "    Coarse grid extent [Xmin, Ymin, Xmax, Ymax]: [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732]\n",
      "    Fine grid extent [Xmin, Ymin, Xmax, Ymax]:   [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732]\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "  Creating CF-netCDF File.\n",
      "    Esri PE String: PROJCS[\"unnamed\",GEOGCS[\"GCS_Sphere\",DATUM[\"D_Sphere\",SPHEROID[\"unnamed\",6370000.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-97.0],PARAMETER[\"Standard_Parallel_1\",30.0],PARAMETER[\"Standard_Parallel_2\",60.0],PARAMETER[\"Latitude_Of_Origin\",41.4710083007812],UNIT[\"Meter\",1.0]]\n",
      "    Map Projection of input raster : lambert_conformal_conic\n",
      "    Starting Process: Building to XMap/YMap\n",
      "    Conversion of input raster to XMap/YMap completed without error.\n",
      "  netCDF global attributes set after  0.01 seconds.\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.13 seconds.\n",
      "  Deriving geocentric coordinates on routing grid from bilinear interpolation of geogrid coordinates.\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "  Creating CF-netCDF File.\n",
      "    Esri PE String: PROJCS[\"unnamed\",GEOGCS[\"GCS_Sphere\",DATUM[\"D_Sphere\",SPHEROID[\"unnamed\",6370000.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-97.0],PARAMETER[\"Standard_Parallel_1\",30.0],PARAMETER[\"Standard_Parallel_2\",60.0],PARAMETER[\"Latitude_Of_Origin\",41.4710083007812],UNIT[\"Meter\",1.0]]\n",
      "    Map Projection of input raster : lambert_conformal_conic\n",
      "    Starting Process: Building to XMap/YMap\n",
      "    Conversion of input raster to XMap/YMap completed without error.\n",
      "    Proceeding to add LATITUDE and LONGITUDE variables after     0.01 seconds.\n",
      "  netCDF global attributes set after  0.01 seconds.\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "    Process: landuse written to output netCDF.\n",
      "Terrain processing step initiated...\n",
      "    Using WhiteboxTools v1.4.0 by Dr. John B. Lindsay (c) 2017-2020\n",
      "    Depression Filling algorithm: Whitebox Fill Depressions.\n",
      "    Process: TOPOGRAPHY written to output netCDF.\n",
      "    Coerced 21 0-value flow direction cells to flow off of the grid.\n",
      "    Process: FLOWDIRECTION written to output netCDF.\n",
      "    Process: FLOWACC written to output netCDF.\n",
      "    Flow accumulation will be thresholded to build channel pixels.\n",
      "    Process: CHANNELGRID written to output netCDF.\n",
      "    Process: STREAMORDER written to output netCDF.\n",
      "    Process: RETDEPRTFAC written to output netCDF.\n",
      "    Process: OVROUGHRTFAC written to output netCDF.\n",
      "    Process: LKSATFAC written to output netCDF.\n",
      "Terrain processing step completed without error in  0.05 seconds.\n",
      "    Forecast points provided and basins being delineated.\n",
      "    Process: frxst_pts written to output netCDF.\n",
      "    Process: basn_msk written to output netCDF.\n",
      "    Channelgrid will not be masked to basins.\n",
      "    Built forecast point outputs in  0.11 seconds.\n",
      "Beginning to build 2D groundwater basin inputs\n",
      "  Building groundwater inputs using FullDom LINKID local basins\n",
      "    Generating LINKID grid for building local sub-basins.\n",
      "Finished building fine-grid groundwater basin grids in  0.02 seconds\n",
      "Beginning to build coarse-grid groundwater basins and parameters\n",
      "    Found 72 basins in the watershed grid\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 1000.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "    Found 69 basins (potentially including nodata values) in the file after resampling to the coarse grid.\n",
      "    GDAL Data type derived from input array: 5 (int32)\n",
      "  Creating CF-netCDF File.\n",
      "    Esri PE String: PROJCS[\"unnamed\",GEOGCS[\"GCS_Sphere\",DATUM[\"D_Sphere\",SPHEROID[\"unnamed\",6370000.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-97.0],PARAMETER[\"Standard_Parallel_1\",30.0],PARAMETER[\"Standard_Parallel_2\",60.0],PARAMETER[\"Latitude_Of_Origin\",41.4710083007812],UNIT[\"Meter\",1.0]]\n",
      "    Map Projection of input raster : lambert_conformal_conic\n",
      "    Starting Process: Building to XMap/YMap\n",
      "    Conversion of input raster to XMap/YMap completed without error.\n",
      "  netCDF global attributes set after  0.01 seconds.\n",
      "    NC dimensions: 16, 15\n",
      "    GWBUCKS array dimensions: 16, 15\n",
      "    Process: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/scratchdir/GWBASINS.nc completed without error\n",
      "    Finished building groundwater grid file in  0.01 seconds\n",
      "    Calculating size and ID parameters for basin polygons.\n",
      "    Created output bucket parameter table (.nc): /home/docker/GIS_Training/Outputs/Gridded_no_lakes/scratchdir/GWBUCKPARM.nc.\n",
      "  Finished building groundwater bucket parameter table in  0.00 seconds.\n",
      "Finished building groundwater parameter files in  0.03 seconds\n",
      "Built output .zip file in  0.49 seconds.\n",
      "Process completed in 0.49 seconds.\n",
      "Archive:  /home/docker/GIS_Training/Outputs/Gridded_no_lakes/Gridded_r4_t20_frxst_mask.zip\n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/GWBASINS.nc  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/GEOGRID_LDASOUT_Spatial_Metadata.nc  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/Fulldom_hires.nc  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/GWBUCKPARM.nc  \n"
     ]
    }
   ],
   "source": [
    "# Create the output directory and ensure it is empty\n",
    "! mkdir -p /home/docker/GIS_Training/Outputs/Gridded_no_lakes\n",
    "! rm -rf /home/docker/GIS_Training/Outputs/Gridded_no_lakes/*\n",
    "\n",
    "# Run the GIS Pre-processing script (with line-breaks)\n",
    "! python Build_Routing_Stack.py \\\n",
    "    -i /home/docker/wrf-hydro-training/example_case/Gridded_no_lakes/DOMAIN/geo_em.d01.nc \\\n",
    "    --CSV /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv \\\n",
    "    -d /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif \\\n",
    "    -R 4 \\\n",
    "    -t 20 \\\n",
    "    -o /home/docker/GIS_Training/Outputs/Gridded_no_lakes/Gridded_r4_t20_frxst_mask.zip\n",
    "\n",
    "# Unzip the directory in-place\n",
    "! unzip /home/docker/GIS_Training/Outputs/Gridded_no_lakes/Gridded_r4_t20_frxst_mask.zip -d /home/docker/GIS_Training/Outputs/Gridded_no_lakes/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### C. Reach-based channel routing configuration with forecast points\n",
    "\n",
    "This domain is a reach-based routing configuration with gages, no masking, no lakes. Regridding factor of 4 and threshold of 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Script initiated at Tue Nov  3 00:15:25 2020\n",
      "  Parameter values that have not been altered from script default values:\n",
      "    Using default basin mask setting: False\n",
      "    Using default OVROUGHRTFAC parameter value: 1.0\n",
      "    Using default RETDEPRTFAC parameter value: 1.0\n",
      "  Values that will be used in building this routing stack:\n",
      "    Input WPS Geogrid file: /home/docker/wrf-hydro-training/example_case/Reach/DOMAIN/geo_em.d01.nc\n",
      "    Forecast Point CSV file: /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv\n",
      "    Mask CHANNELGRID variable to forecast basins?: False\n",
      "    Create reach-based routing (RouteLink) files?: True\n",
      "    Lake polygon feature class: \n",
      "    Input high-resolution DEM: /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif\n",
      "    Regridding factor: 4\n",
      "    Stream initiation threshold: 20\n",
      "    OVROUGHRTFAC parameter value: 1.0\n",
      "    RETDEPRTFAC parameter value: 1.0\n",
      "    Input channel initiation start point feature class: None\n",
      "    Input groundwater basin polygons: None\n",
      "    Output ZIP file: /home/docker/GIS_Training/Outputs/Reach/Reach_r4_t20_frxst.zip\n",
      "  Running Process GEOGRID function\n",
      "    Forecast points provided.\n",
      "  Reach-based routing files will be created.\n",
      "WPS netCDF projection identification initiated...\n",
      "    Map Projection: Lambert Conformal Conic\n",
      "    Using MOAD_CEN_LAT for latitude of origin.\n",
      "    Using Standard Parallel 2 in Lambert Conformal Conic map projection.\n",
      "    Geo-referencing step completed without error in  0.04 seconds.\n",
      "  Building sub-grid of model grid.\n",
      "    Original grid spacing dx=1000.0, dy=-1000.0\n",
      "    Original grid size: rows=16, cols=15\n",
      "    New grid spacing: dx=250.0, dy=-250.0\n",
      "    New dimensions: rows=64, cols=60\n",
      "    Created projection definition from input NetCDF GEOGRID file.\n",
      "    Proj4: +proj=lcc +lat_0=41.4710083007812 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs\n",
      "    Coarse grid GeoTransform: 1841999.0194101864 1000.0 0 278495.814356732 0 -1000.0\n",
      "    Coarse grid extent [Xmin, Ymin, Xmax, Ymax]: [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732]\n",
      "    Fine grid extent [Xmin, Ymin, Xmax, Ymax]:   [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732]\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "  Creating CF-netCDF File.\n",
      "    Esri PE String: PROJCS[\"unnamed\",GEOGCS[\"GCS_Sphere\",DATUM[\"D_Sphere\",SPHEROID[\"unnamed\",6370000.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-97.0],PARAMETER[\"Standard_Parallel_1\",30.0],PARAMETER[\"Standard_Parallel_2\",60.0],PARAMETER[\"Latitude_Of_Origin\",41.4710083007812],UNIT[\"Meter\",1.0]]\n",
      "    Map Projection of input raster : lambert_conformal_conic\n",
      "    Starting Process: Building to XMap/YMap\n",
      "    Conversion of input raster to XMap/YMap completed without error.\n",
      "  netCDF global attributes set after  0.01 seconds.\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.13 seconds.\n",
      "  Deriving geocentric coordinates on routing grid from bilinear interpolation of geogrid coordinates.\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "    GDAL Data type derived from input array: 6 (float32)\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "  Creating CF-netCDF File.\n",
      "    Esri PE String: PROJCS[\"unnamed\",GEOGCS[\"GCS_Sphere\",DATUM[\"D_Sphere\",SPHEROID[\"unnamed\",6370000.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-97.0],PARAMETER[\"Standard_Parallel_1\",30.0],PARAMETER[\"Standard_Parallel_2\",60.0],PARAMETER[\"Latitude_Of_Origin\",41.4710083007812],UNIT[\"Meter\",1.0]]\n",
      "    Map Projection of input raster : lambert_conformal_conic\n",
      "    Starting Process: Building to XMap/YMap\n",
      "    Conversion of input raster to XMap/YMap completed without error.\n",
      "    Proceeding to add LATITUDE and LONGITUDE variables after     0.01 seconds.\n",
      "  netCDF global attributes set after  0.01 seconds.\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 250.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "    Process: landuse written to output netCDF.\n",
      "Terrain processing step initiated...\n",
      "    Using WhiteboxTools v1.4.0 by Dr. John B. Lindsay (c) 2017-2020\n",
      "    Depression Filling algorithm: Whitebox Fill Depressions.\n",
      "    Process: TOPOGRAPHY written to output netCDF.\n",
      "    Coerced 21 0-value flow direction cells to flow off of the grid.\n",
      "    Process: FLOWDIRECTION written to output netCDF.\n",
      "    Process: FLOWACC written to output netCDF.\n",
      "    Flow accumulation will be thresholded to build channel pixels.\n",
      "    Process: CHANNELGRID written to output netCDF.\n",
      "    Process: STREAMORDER written to output netCDF.\n",
      "    Process: RETDEPRTFAC written to output netCDF.\n",
      "    Process: OVROUGHRTFAC written to output netCDF.\n",
      "    Process: LKSATFAC written to output netCDF.\n",
      "Terrain processing step completed without error in  0.05 seconds.\n",
      "    Forecast points provided and basins being delineated.\n",
      "    Process: frxst_pts written to output netCDF.\n",
      "    Process: basn_msk written to output netCDF.\n",
      "    Channelgrid will not be masked to basins.\n",
      "    Built forecast point outputs in  0.11 seconds.\n",
      "    Routing table will be created...\n",
      "        Stream to features step complete.\n",
      "        Warning: Number of unique IDs exceeds limit of 16-bit unsigned integer type. Not all reaches may be converted to stream vectors. Check output carefully.\n",
      "        Found 69 unique IDs in stream vector layer.\n",
      "        Eliminating 3 IDs in LINKID grid that could not be resolved in stream vector layer.\n",
      "          [43, 44, 65]\n",
      "        Adding forecast points:LINKID association.\n",
      "        Found 3 forecast point:LINKID associations.\n",
      "  All dictionaries have been created in  0.00 seconds.\n",
      "  Fields have been added to the shapefile.\n",
      "    Projection defined for input vector layer in 0.00 seconds.\n",
      "        Starting to fill in routing table NC file.\n",
      "        Done writing NC file to disk.\n",
      "    Routing table created without error.\n",
      "Reach-based routing inputs generated in 0.07 seconds.\n",
      "Beginning to build 2D groundwater basin inputs\n",
      "  Building groundwater inputs using FullDom LINKID local basins\n",
      "    Generating LINKID grid for building local sub-basins.\n",
      "Finished building fine-grid groundwater basin grids in  0.01 seconds\n",
      "Beginning to build coarse-grid groundwater basins and parameters\n",
      "    Found 72 basins in the watershed grid\n",
      "    Raster resampling initiated...\n",
      "    The High-resolution dataset will be 1000.0m\n",
      "    Projected input raster to model grid in  0.02 seconds.\n",
      "    Found 69 basins (potentially including nodata values) in the file after resampling to the coarse grid.\n",
      "    GDAL Data type derived from input array: 5 (int32)\n",
      "  Creating CF-netCDF File.\n",
      "    Esri PE String: PROJCS[\"unnamed\",GEOGCS[\"GCS_Sphere\",DATUM[\"D_Sphere\",SPHEROID[\"unnamed\",6370000.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",-97.0],PARAMETER[\"Standard_Parallel_1\",30.0],PARAMETER[\"Standard_Parallel_2\",60.0],PARAMETER[\"Latitude_Of_Origin\",41.4710083007812],UNIT[\"Meter\",1.0]]\n",
      "    Map Projection of input raster : lambert_conformal_conic\n",
      "    Starting Process: Building to XMap/YMap\n",
      "    Conversion of input raster to XMap/YMap completed without error.\n",
      "  netCDF global attributes set after  0.01 seconds.\n",
      "    NC dimensions: 16, 15\n",
      "    GWBUCKS array dimensions: 16, 15\n",
      "    Process: /home/docker/GIS_Training/Outputs/Reach/scratchdir/GWBASINS.nc completed without error\n",
      "    Finished building groundwater grid file in  0.01 seconds\n",
      "    Calculating size and ID parameters for basin polygons.\n",
      "    Created output bucket parameter table (.nc): /home/docker/GIS_Training/Outputs/Reach/scratchdir/GWBUCKPARM.nc.\n",
      "  Finished building groundwater bucket parameter table in  0.00 seconds.\n",
      "Finished building groundwater parameter files in  0.03 seconds\n",
      "Built output .zip file in  0.55 seconds.\n",
      "Process completed in 0.55 seconds.\n",
      "Archive:  /home/docker/GIS_Training/Outputs/Reach/Reach_r4_t20_frxst.zip\n",
      "  inflating: /home/docker/GIS_Training/Outputs/Reach/streams.prj  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Reach/streams.shp  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Reach/streams.dbf  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Reach/GWBASINS.nc  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Reach/streams.shx  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Reach/GEOGRID_LDASOUT_Spatial_Metadata.nc  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Reach/Fulldom_hires.nc  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Reach/GWBUCKPARM.nc  \n",
      "  inflating: /home/docker/GIS_Training/Outputs/Reach/Route_Link.nc  \n"
     ]
    }
   ],
   "source": [
    "# Create the output directory and ensure it is empty\n",
    "! mkdir -p /home/docker/GIS_Training/Outputs/Reach\n",
    "! rm -rf /home/docker/GIS_Training/Outputs/Reach/*\n",
    "\n",
    "# Run the GIS Pre-processing script (with line-breaks)\n",
    "! python Build_Routing_Stack.py \\\n",
    "    -i /home/docker/wrf-hydro-training/example_case/Reach/DOMAIN/geo_em.d01.nc \\\n",
    "    -r True \\\n",
    "    --CSV /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv \\\n",
    "    -d /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif \\\n",
    "    -R 4 \\\n",
    "    -t 20 \\\n",
    "    -o /home/docker/GIS_Training/Outputs/Reach/Reach_r4_t20_frxst.zip\n",
    "\n",
    "# Unzip the directory in-place\n",
    "! unzip /home/docker/GIS_Training/Outputs/Reach/Reach_r4_t20_frxst.zip -d /home/docker/GIS_Training/Outputs/Reach/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### You have reached the end of Lesson S2"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
